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In the quasi-stationary states of the Hamiltonian Mean-Field model, we numerically compute 
correlation functions of momenta and diffusion of angles with homogeneous initial conditions. This 
is an example, in a N-body Hamiltonian system, of anomalous transport properties characterized by 
non exponential relaxations and long-range temporal correlations. Kinetic theory predicts a striking 
transition between weak anomalous diffusion and strong anomalous diffusion. The numerical results 
are in excellent agreement with the quantitative predictions of the anomalous transport exponents. 
Noteworthy, also at statistical equilibrium, the system exhibits long-range temporal correlations: 
the correlation function is inversely proportional to time with a logarithmic correction instead of 
the usually expected exponential decay, leading to weak anomalous transport properties. 

PACS numbers: 

05. 45. -a Nonlinear dynamics and nonlinear dynamical systems, 
05.70.Ln Noncquilibrium and irreversible thermodynamics. 



I. INTRODUCTION 

Recently, a new light was shed on long-range interacting systems The first reason is that a mathematical 
characterization Q and the study of several simple models have completely clarified the inequivalence of ensembles 
that might exists between the microcanonical and the canonical ensembles 0, H| . The second is the appearance of 
a very useful technique, namely the large deviation theory, to compute the microcanonical number of microstates 
and thus the associated microcanonical entropy Q . The third is a classification of all possible situations of ensemble 
inequivalence Q. The last, but not the least, reason is the understanding that the broad spectrum of applications (self- 
gravitating Q and Coulomb systems, vortices in two-dimensional fluid mechanics, wave-particles interaction, trapped 
charged particles, ...) [l[ should be considered simultaneously since significant advances were performed independently 
in the different domains. However as usual in Physics, the study of simple models is of particular interest not only for 
pedagogical properties, but also for testing ideas that might be derived analytically and verified numerically without 
very expensive simulations. 

We consider here the Hamiltonian Mean Field (HMF) model, which is considered as the paradigmatic dynamical 
model for long-range interacting systems. This model [1, [t| [l(| El consists of N particles moving on the unit circle, 
and is described by the Hamiltonian 

N N N 
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where 8j is the angle of j-th particle and pj its conjugate momentum. Using a change of the time unit, the prefactor 
l/N of the second term is added in order to get an extensive energy Q. Thus, in the limit N — > oo, the appropriate 
mean-field scaling is obtained for the statistical mechanics. Studies of the HMF model have been recently reinforced 
by the discovery of its link with the Colson-Bonifacio model for the single-pass free electron laser Q . 

Within this model, a striking disagreement was reported between the canonical statistical mechanics predictions 
and time averages of constant energy molecular dynamics simulations [n], IH|. As the model has only a second 
order phase transition [f| at the critical energy density U c = 3/4, the possibility that the origin might lead to an 
inequivalence between canonical and microcanonical statistical mechanics can be excluded j6[. Moreover, recently, it 
has been shown unambiguously that the microcanonical entropy leads to the same predictions than the canonical free 
energy [B[ . Very interesting results about the behavior of such a system in contact with a thermal bath has however 
been recently reported [lj, LL4[ . 

The origin of the apparent disagreement comes from a particularly slow dynamical evolution of this long-range 
system. Indeed, in Hamiltonian systems with long-range interactions, systems are sometimes trapped in quasi- 
stationary states (QSS) before going to equilibrium. Examples of such QSS were found in a 1-dimcnsional self- 
gravitating system [15[ and in the HMF model [l6[. The trapping time diverges algebraically in the limit N — > oo 
and, hence, time averages disagree with canonical averages if the computing time is not long enough. 
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To understand th e dy namics during such a long period, QSS were interpreted as stable stationary states of the 
Vlasov equation [13, [Hi that can be derived from the Hamiltonian dynamics. The Vlasov equation, which governs 1- 
particle distribution function is indeed exact (l9l | in the limit N — > oo, but only approximate for a finite system: finite 
size effects drive indeed the system from the Vlasov stable stationary state to the Boltzmann equilibrium. Recently, 
Caglioti and Rousset [2(| proved for a wide class of potentials which includes the HMF case, that iV particles starting 
close to a Vlasov stable stationary state remain close to it during a time scale proportional at least to V 1 / 8 . The 
result is consistent with numerical results which state that the lifetime of QSS scales like N 1 - 7 fl(| . 

Using a kinetic approach which goes beyond the above Vlasov interpretation, the correlation function of momenta 
was recently derived [2l|, [22j with the following assumptions: (i) a finite but large enough number of particles, (ii) a 
homogeneous distribution of angles, and (hi) a system in a (quasi-)stationary state. As shown in Refs. [IE EH, the 
latter condition amounts to consider initial distributions of momenta fo(p) satisfying the inequality 

1 ar«M dP>0 . ( 2 ) 

2 J-oo P 

This condition has been derived for the linear [23| and formal [13] stability of the distribution /o (see also Ref. [24| 
for another derivation). This condition defines a critical energy U* which is, in general, different from the critical 
energy U c = 3/4 where the second order phase transition is located. However, as expected, both values coincide for 
a gaussian distribution fo(p)- Above theory is expected to be valid in the time interval 1 <C r -C N, where r = t/N 
is the appropriate rescaled time. 

Among the main predictions resumed in Table HI one might emphasize that distributions /o(p) with algebraic tails 
were proved to have a correlation function of momenta C p (r) with an algebraic decay in the long-time regime. Striking 
algebraic large time behaviors for momentum autocorrelations had been first numerically observed in Refs. [H|,[25j]. On 
the contrary, distributions with stretched exponential tails correspond to correlation functions inversely proportional 
to time with a logarithmic correction. It is also important to stress that gaussian distributions, which corresponds 
to 6 = 2 in the stretched exponential case, leads to a long-time correlation of In r/r instead of the usual exponential 
decay in the stable, supercritical energy regime U > U* = U c , although the initial distribution is at equilibrium. 
The origin of this unusual long-time momentum correlations does not depend on the center part of the momentum 
distributions fo(p) but on its tails. One might understand physically this behavior, by the fact that particles located 
in these tails move almost freely, and hence yield long-time correlations. 
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TABLE I: Asymptotic forms of initial distributions fo(p), and theoretical predictions of correlation functions C p (t) and the 
diffusion <jg (r) in the long-time regime. Asymptotic forms of the distribution and the predictions are assumed and predicted in 
the limits \p\ —> oo and r — > oo respectively, where r = t/N is a rescaled time. The exponent a is given as ot = (u — 3)/(u + 2) . 
See Ref. [23] for details. 

In these (quasi-)stationary states, the theoretical law for the diffusion of angles cr|(r) has been also derived. The 
predictions [22, [26| for the diffusion properties are listed in Table HI They clarify the highly debated disagreement 
between different numerical simulations reporting either anomalous [23 ] or normal [l6j diffusion, in particular by 
delineating the time regime for which such anomalous behavior should occur. We briefly recall that when the moment 
of order n of the distribution scales like r' 1 / 2 at large time, such a transport is called normal. However, anomalous 
transport pll. l29l. [30I] . where moments do not scale as in the diffusive case, were reported in some stochastic models, in 
continuous time random walks (Levy walks) , and for systems with a lack of stationarity of the corresponding stochastic 
process [3l|. When the distribution fo(p) is changed within the HMF model, a transition between weak anomalous 
diffusion (normal diffusion with logarithmic corrections) and strong anomalous diffusion is thus predicted. From the 
physical point of view, as particles with large momentum p fly very fast in comparison to the typical time scales of 
the fluctuations of the potential, they are subjected to a very weak diffusion and thus maintain their large momentum 
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during a very long time. A thick distribution of waiting time with a large momentum explains the anomalous diffusion. 
From a mathematical point of view, these behaviors are linked to the non exponential relaxation of the Fokker-Planck 
equation describing the diffusion of momenta, leading to long-range temporal correlations (22J. This mechanism is 
new in the context of kinetic theory. However, similar Fokker-Planck equations, with a rapidly vanishing diffusion 
coefficients obtained by other physical mechanisms, have been studied in several frameworks [32l. l33l. |34| . 

The first purpose of this article is to numerically check the theoretical predicted correlation functions for power-tail 
and gaussian distributions by using accurate numerical simulations. The other is to clarify whether diffusion is normal 
or anomalous, which depends on the choice of the initial distribution fo(p)- 

The article is organized as follows. Some useful quantities are introduced in Section [TTJ In Sections [HT] and IIV1 
we respectively focus on initial distributions with power-law and gaussian tails. In each section, we first check the 
stationarity and the stability following the method developed in Ref. [f7| and determine the time region of the QSS. We 
also study carefully the correlation function and the diffusion comparing them with theoretical predictions. Finally, 
section M concludes the discussion. 



II. QUANTITIES OF INTEREST AND NUMERICAL PROTOCOL 

In order to check the stationarity and the stability of an initial distribution / (p) , we study the temporal evolutions 
of several macrovariables: 

• the magnetization defined as the modulus M of the vector M = (M x , M y ), where both components are defined 
as M x = (cos9) N and M y = (s'md) N . The bracket (-) N represents the average over all particles, for instance 
(cos9) N = (J^Li cosdj)/N. Note that the magnetization M is constant if the system is stable stationary. 

• the moments of the 1-body distribution function f(9,p 7 t). As explained in details in Ref. [l7| ] . the stationarity 
of the 1-body distribution f(9,p,t) implies the stationarity of the individual energy distribution / e (e,i), where 
e = p 2 /2 — M x cos9 — M y sm9. Moreover, the stationarity of f e (e,i) implies the stationarity of all moments 
l^n = (e")jv- As the stationarity of the moment is a necessary condition for stability, vanishing derivatives 
fi n = d/z ra /di, for n = 1,2 and 3, would suggest that the system is in a (quasi-)stationary state, while large 
derivatives clearly indicate a non-stationary state. In addition, the stability is suggested if a state stays stationary 
for a long period. 

While checking the stationarity and the stability, we identify a time region where the system is in the QSS, during 
which we observe the correlation function of momenta C p (r) = (p(r)p(0)) N and the diffusion of angles crg(r) = 
([9(t) — 9(0)} 2 ) N . The latter quantity can be rewritten as follows 



N 2 



u 



dn / dr 2 (p(n>(r 2 )} w =2 da dr 2 (p(a + t 2 ) P (t 2 )) n , (3) 



where the factor 1/N 2 comes from the time rescaling r = t/N, while the new variable s = t\ — r 2 was introduced to 
take advantage of the division of the square domain into two isoscale triangles corresponding to s > and s < 0. In 
the (quasi-)stationary states, the integrand (p(s + T 2 )jo(T 2 )) at does not depend on r 2 (the QSS evolve on a time scale 
much larger than N) and hence diffusion can be simplified jl6| by using the correlation function as 



N 2 



(r-s)C p (s) ds. (4) 



We numerically performed the temporal evolution of the canonical equations of motion by using a 4-th order 
symplectic integrator [35|, HH with a time step At = 0.2 and a total momentum set to zero. Initial values of angles 
being randomly chosen from a homogeneous distribution, the magnetization M is hence of order l/yiV. Omitting 
this vanishing value of M, the energy density U = K + (1 — M 2 )/2 where K is the kinetic energy density can thus 
be well approximated by the kinetic energy density K as U = K + 1/2. To characterize the simulations, the only 
remaining point is the initial distribution of momenta: in the following sections, as anticipated, we will carefully study 
distributions with power-law and gaussian tails. 
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III. POWER-LAW TAILS 



A. Initial distribution 



In this section, we consider the initial distribution 

/o(p) = 



i + Ip/poI"' 



(5) 



whose power-law tails are characterized by the exponent v. The unity, added in the denominator to avoid the 
divergence at the origin p = 0, does not affect neither the asymptotic form, nor the theoretical predictions. The 

1 /2 

parameter po is directly determined by the kinetic energy density as po = (2_ftTsin(37r/j/)/sin(7r/^)) ' , while the 
normalization factor is 



v ( sin 3 (7r/f) 



2tt \2KsSi\{Zit/v) 



1/2 



(0) 



From the stability criterion one gets that this initial state is Vlasov stable when the kinetic energy density satisfies 
the condition 



K > 



1 sin(7r/^) 
4 sin(37r/z/) 



(7) 



One thus gets a dynamical critical energy U* = 0.75, 0.625 and 0.60355 ... for v = 4, 6 and 8 respectively. In the 
rest of this section, we set the exponent v to 8. 

B. Stationarity and stability checks 

Let us numerically check the stationarity and the stability of these states; in particular, it will clarify the time 
region of existence of the QSS. Figure Q] presents the temporal evolution of M and (i n (n = 1,2,3) for the unstable 
(U = 0.6 < U*) and stable (U = 0.7 > U*) cases. In both cases, the magnetization M eventually goes toward the 
equilibrium value M eq , indicated by horizontal lines. The three quantities fi n have vanishing small fluctuations around 
zero except during the time interval 0.0005 < r < 0.003 for the unstable case. In the unstable case, the system is first 
in an unstable stationary state (r < 0.0005), before becoming non-stationary (0.0005 < r < 0.003) and finally reaches 
stable stationary states (t > 0.003). On the other hand, in the stable case, the stable stationarity holds throughout 
the computed time. 
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FIG. 1: Stationarity check for initial distributions with power- law tails. Note the logarithmic scale for the rescaled time 
r = t/N. Panel (a) presents the unstable case U = 0.6 while panel (b) the stable one U = 0.7. The three curves fin (n — 1, 2, 3) 
are reported in both panels. Their vertical scales are multiplied by 100 for graphical purposes. Curves and horizontal lines 
indicated by symbols M and M eq represent respectively temporal evolutions of the magnetization and its equilibrium value. 
All numerical curves are obtained by averaging 20 different numerical simulations for N = 10 4 . 



In the stable case, the magnetization M stays around zero before taking off around r = 20 to reach the equilibrium 
value M eq . The fluctuation level of fi n increases around the take-off time r = 20, but the increase does not imply 
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any non-stationarity of the system, since the fluctuation level is 10 times smaller than the corresponding one in the 
non-stationary time region of the unstable case. The nonzero magnetization might be at the origin of the larger 
fluctuations than in the zero magnetization cases, since the former has a phase and an individual energy e which 
depends not only on the modulus M but also on the phase. 



C. Check of the theoretical prediction 

In the stable case (U = 0.7), we perform numerical computations for N = 10 3 , 10 4 , 2.10 4 and 5.10 4 , and average 
over 20, 20, 10 and 5 sample orbits respectively. Temporal evolutions of magnetization M are shown in Fig. [DJa), and 
M takes off toward the equilibrium value M eq around r 2 = 1, 20 and 50 for N = 10 3 , 10 4 and 2.10 4 respectively. The 
take-off time defines the end of applicable time region of the theory since the homogeneous assumption (ii) breaks. 
Note that no take-off time appears in the case N = 5.10 4 , within the computed time interval. 

The theory predicts (see Table Q] for v = 8) that the correlation function decays algebraically with the exponent 
— 1/2, i.e. C p (t) ~ t" 1 / 2 , up to the take-off time r 2 . According to Fig. EJb), the theoretical prediction agrees well 
with numerical computations in the intermediate time region t\ < t < T2, where n = 2 for any value N. This is 
expected since, on the one hand, the short-time region r < t\ is out of the time domain of application since the theory 
gives asymptotic estimates. The time t\ is marked as a long vertical line in Fig. [2jb) to clearly indicate the start of 
the applicable time domain. Although the quantity t s ~ 0.005 is not derived theoretically, the straight lines with the 
slope —1/2 in Fig. [^b), representing (t/t s ) -1 / 2 , emphasizes the agreement of the predicted exponent. 

Introducing the expression of the correlation function in relation leads to the law <Jq(t) ~ t 3 / 2 : Figure [2Jc), in 
which the four curves for the four different values of N almost collapse, attests also the validity of this prediction in 
the intermediate time region t\ < t < t 2 . It is possible to confirm more precisely that the diffusion exponent is 3/2 
by introducing the instantaneous exponent 7 [371 ] defined as 

djna|(r) 1 da 2 e (r) 
' dlnr ct 2 (t) dlnr ' 1 ' 

The instantaneous exponent 7, shown in Fig.[2jd), goes down and once crosses 3/2. However, 7 comes back and stays 
around 3/2 in the time interval t\ < t < r 2 . Above result confirms therefore unambiguously that the diffusion is 
anomalous, namely superdiffusive, in the intermediate QSS time interval as predicted by the theory [22j|. 

The temporal evolution of 7 was also recently discussed by Antoniazzi et al. (38j . and was shown to monotonically 
decrease toward 1. The difference has two different origins: First, Antoniazzi et al considered non-homogeneous 
initial distribution of angles, which are out of the applicable range of the theory tested here. Second they considered 
a waterbag initial distribution of momenta, which does not have tails initially, although tails develop of course as soon 
as the time is slightly positive. As the theory states that the asymptotic law for diffusion is determined by the tails of 
the initial distribution of momenta, there is no contradiction that the temporal evolution of 7 is different. A similar 
remark applied with the out-of-equilibrium initial distribution discussed by Moyano and Anteneodo (37j . 

For the power-tail initial distributions, the theoretical predictions are essentially good, but not exact. We first note 
that the increase of N does not affect neither the correlation function, nor the diffusion, at least for 10 4 < N < 5.10 4 
(the case N = 10 3 has been excluded since no validity time region t\ < r < t 2 appears). In the numerical results, the 
slope of the diffusion 7 is not 1.5 but belongs to [1.44, 1.48]. The relative discrepancy is thus at most of 4 percents. 
There are two possibilities to understand this small discrepancy: (a) the lack of the samples, or (b) the lack of 
stationarity which is the assumption (iii) of the theory. We will discuss on the origin of these discrepancies at the end 
of the next section. 



IV. GAUSSIAN DISTRIBUTION 



A. Initial distribution 



In this section, we consider the gaussian initial distribution 

where the initial temperature T is determined from the energy density as T = 2K = 2U — 1. The dynamical critical 
energy of this gaussian distribution coincides with the critical energy of the second order phase transition U c = 3/4. As 
the distribution of angles is homogeneous, the system is therefore at equilibrium for any supercritical energy U > U c . 
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FIG. 2: Check of the theoretical prediction for stable initial distributions with power-law tails, in the case U = 0.7. Points are 
numerically obtained by averaging 20, 20, 10 and 5 realizations for N = 10 3 , 10 4 , 2.10 4 and 5.10 4 respectively. Panel (a) shows 
temporal evolution of magnetization in the scale time r = t/N. The take-off times of M are estimated as T2 = 1,20 and 50 
for N = 10 3 , 10 4 and 2.10 4 respectively, and T2 are marked in panel (b) and (d). No take-off for N = 5.10 4 is observed in this 
computing time. The horizontal line represents the equilibrium value of M. In panel (b), four curves represent the correlation 
functions of momenta, while the straight lines with the slope —1/2 represent the theoretical prediction. The curves and the 
lines are multiplied from the original vertical values by 2,4 and 8 for N = 10 4 ,2.10 4 and 5.10 4 for graphical purposes. The 
vertical line indicates n = 2 from which the valid time region of the theory starts. Similarly, panel (c) presents the diffusion 
of angles, while the straight line with the slope 3/2 is theoretically predicted. The four curves for the four different values of 
N are reported and almost collapse. Finally, panel (d) shows the temporal evolution of the instantaneous exponent 7, defined 
in Eq. ([8]), and 7 stays around the theoretically predicted value 3/2 in the time region n < r < r2. The values 0.1, 0.2 and 0.3 
are added in vertical values for N = 10 4 ,2.10 4 and 5.10 4 respectively for graphical purposes. 



B. Stationarity and stability checks 

The stationarity and stability arc checked as in Sec. IIII Bl by considering temporal evolutions of magnetization and 
the derivatives of moments fi n shown in Fig. [3] The scenario of relaxation of this initial distribution with power-law 
tails is very similar. In the unstable case (JJ = 0.7 < U c ), the system reaches stable stationary states after experiencing 
unstable stationary and non-stationary states. In the stable case (U = 0.8 > U c ), the state is stable stationary in the 
whole time domain since it is initially at equilibrium. 



C. Check of the theoretical prediction 

Let us focus on the stable case U = 0.8 with N = 10 4 . The correlation function obtained numerically, and shown 
in Fig. [3Ja), is in good agreement with the theoretical prediction (In r)/r in the long-time region r > t\ — 1 if we 
accept the second scaling of time as r — > t/t s with r s = 0.2. As already mentioned in Sec. IIII Cl the second scaling is 
not provided by the theory, while the asymptotic theoretical estimate is out of applicability in the short time domain 
r < T\ . We would like also to stress that the logarithmic correction makes the prediction more precise rather than a 
simple algebraic decay 1/r. 
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FIG. 3: Stationarity check for gaussian initial distributions with iV = 10 4 . Note the logarithmic scale for the rescaled time 
r = t/N. Panel (a) presents the unstable case U = 0.7 while panel (b) the stable one U = 0.8. The three curves fi n (n — 1, 2, 3) 
are reported in both panels. Their vertical scales are multiplied by 100 for graphical purposes. Curves indicated by symbol 
M represent the temporal evolutions of the magnetization. In panel (a), the horizontal line indicated by M eq represents the 
equilibrium value, while in panel (b), the equilibrium value is zero. All numerical curves are obtained by averaging 20 different 
numerical simulations. 



The correlation function can thus be approximated as 

C p (0) 



C p (t) = 



— In — 



if r < T\ 

if T > T\ 



(10) 



where the short time value has to be C p (0) = (p 2 (0)) Ar = 2K = 0.6, while c = 0.85 is obtained by a fitting procedure. 
This approximation of the correlation function and the relation ((?]) leads to the following expression for the diffusion 

( C p (0)t 2 , ifr<n 



-JfT = \ 2C P (0) 



rir ~T 



-2ct s 



In I- 



ln^ 



Inl-lj-rWln^-l 



(11) 



if T > Tl 



Figure [31(b) presents the diffusion obtained numerically. The two straight lines indicating the short- and long-time 
regions shows a very good agreement. The diffusion seems anomalous with an exponent 1.35 in the long-time region. 
Similarly, the instantaneous exponent 7 seems to converge toward 1.35 as shown by Fig. |UJd). However, these 
observations are not accurate, and only due to a long transient, induced by the logarithmic correction. Diffusion is 
essentially proportional to the time r, and hence must be normal in the asymptotic time region. Expression (|11[) 
provides the asymptotic form of the instantaneous exponent 
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ln(r/r s ) 



(12) 



This prediction is in good agreement with numerical results as attested by Fig. |4jd). I n the limit of r — > 00, the 
exponent 7 goes logarithmically toward unity, and we therefore conclude that diffusion is normal although a long 
transient time is necessary to observe it. This is an excellent illustration of the difficulty to get reliable numerical 
estimates for the diffusion exponent. Such a case explains a posteriori the reason of previous disagreement 

nana. 

As predicted by Table QJ the logarithmic correction of the correlation function yields a logarithmic correction of the 
diffusion, so that its asymptotic form should be a"g(T)/N 2 ~ r(lnr) 2 . Figure HJc) confirms this prediction by plotting 
y/ Vg(r) / '(tN 2 ) as a function of lnr: one gets a linear behavior in the long time region r > 1. We thus have confirmed 
the existence of weak anomalous diffusion, i.e. normal diffusion with logarithmic corrections. 

Let us return to the origin of discrepancies for a and 7, discussed at the end of the previous section for the power- 
law tails. It seems natural to exclude the possibility (a), lack of samples, since the same number of orbits, 20, has 
been used in the case N — 10 4 , for both the power-tails and the gaussians, while the latter case agrees extremely 
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FIG. 4: Check of the theoretical prediction for gaussian stable initial distributions in the case U = 0.8 with N = 10 4 . Points 
are numerically obtained by averaging 20 realizations. In panel (a), the symbols show the correlation function of momenta. The 
theoretical prediction, lnr/r, is a better approximation than the simpler law 1/r. Panel (b) presents the diffusion of angles. 
Although the diffusion is normal, the straight line r 1,35 wrongly suggests that it is not. See text for explanations and details. 
Panel (c) shows the quantity \/o~g(T)/(rN 2 ) as a function of lnr to confirm the logarithmic correction of the diffusion. The 
straight line is a guide for the eyes. Finally, panel (d) presents the temporal evolution of the instantaneous exponent 7. The 
dashed line corresponds to relation (|12[1 , 



well with the theoretical predictions, even including the logarithmic correction. This excellent agreement comes from 
the absence of any breaking of theoretical assumptions, since the state is at equilibrium and stationary accordingly. 
Consequently, we can consider that the possibility (b), lack of stationarity, explains the discrepancies of a and 7 for 
power-law tails. 



V. SUMMARY 

We have numerically confirmed the theoretical predictions proposed in Ref. [22j for initial distributions with power- 
law or gaussian tails: correlation function and diffusion are in good agreement with numerical results. Diffusion 
is indeed anomalous superdiffusion in the case of power-law tails, while normal when gaussian. In the latter case, 
the system is at equilibrium, but the diffusion exponent shows a logarithmically slow convergence to unity due to 
a logarithmic correction of the correlation function. This long transient time to observe normal diffusion, even for 
gaussian distribution and at equilibrium, suggests that one should be very careful to decide whether diffusion is 
anomalous or not [3^, |jf| HH ■ 

For the power-law tails initial distribution, the numerically obtained exponent of diffusion is slightly different 
from the theoretical prediction (few percents). As discussed above, this discrepancy comes from the breaking of the 
stationary assumption. The state is only approximately stationary, explaining that the theoretical predictions are not 
exact but only approximate. We stress that in the limit of large N, these states become stationary because their living 
times diverge much faster than N. For the gaussian initial distribution, the state is in equilibrium from the start, and 
stationary even with finite N: hence the theoretical predictions agree extremely well with numerical results. 

In addition, above numerical computations clarify two new points: (i) the time region where the theory is applicable, 
(ii) the second time scaling to fit the correlation function and the diffusion. Both might depend on the degrees of 
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freedom, but the latter, (ii), appears to be not the case for the power-law tails. Obtaining the dependence for the 
gaussian is a future work. 

Finally, let us remark that the scenario of the relaxation described in Refs. [13, EH is confirmed even for initial 
distributions with power-law tails: this had never been tested previously. The scenario claims that the system with 
long-range interactions experiences first a violent relaxation, before the so-called collisional relaxation which drives 
the system toward Boltzmann's equilibrium. In the simulations reported here, non-stationary and stable stationary 
states correspond respectively to the violent and the collisional relaxations. One might also remark that distributions 
with power-law tails induce quasi-stationary states above the dynamical critical energy, while being not a member of 
^-distributions [42|. The latter might be a sufficient condition of QSS, but is definitely not a necessary condition. To 
conclude let us remark that if the results discussed here concerns the simple HMF model, let us mention that it is 
believed to be general for long-range interacting systems [H, [44| . 
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